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1 Introduction 

In this report, a new algorithm for the solution of the steady Euler equations in two- 
dimensional flow is presented. The development of the scheme was motivated by several 
considerations. First, the geometric generality of an unstructured grid was desired to 
allow the ease of treating complex geometries. Second, this generality should not be at 
the expense of accuracy; it was required that the scheme be second-order accurate in 
the steady-state Third, it. was desired to keep the data structure as simple as possible 
to avoid excess overhead in terms of storage, and to minimize communications costs 
on any massively parallel computer that may be used in the future. This requires 
a compact difference stencil, so that the grid connectivity information necessary to 
update the solution at a node is minimal Finally, the use of upwind differencing was 
desired to provide sharp resolution of shocks and to minimize the need for adjusting 
artificial viscosity coefficients, in contrast to some well-established central-difference and 
finite-element schemes. 

The scheme that has emerged from the above requirements is called the Upwind 
Control Volume (UCV) scheme. The ideas for this algorithm have come from two 
independent lines of thought. The first comes from the version of the Lax-^Vendroff 
algorithm developed by Ni (1]. Ni’s interpretation of his algorithm is that it is a form 
of* upwind scheme, in which the flux residual is evaluated on a cell and then distributed 
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preferentially towards the downstream nodes. In one dimension, a scalar Lax- YVendrofT 
algorithm with a Courant number of 1.0 becomes an exact upwind scheme, i.e. it pro- 
duces the exact convection of the initial conditions. However, at smaller values of the 
Courant number the upwind biasing is lessened because it is a second order term. This 
is particularly important in solving systems of equations, such as the Euler equations. 
At a shock, if the Courant number associated with the Ufa characteristic is 1.0, then 
the Courant number of the u a characteristic is very much smaller, and therefore is 
only weakly upwinded. However, the latter characteristic is the one that defines the 
shock. This is the reason that artificial viscosity terms must be added to Ni’s Lax- 
Wendroff calculations to prevent unacceptable oscillations a' shocks. The second line 
of thought comes from the ideas of Moore [2j, who has developed a three-dimensional 
elliptic pressure-correction Navier-Stokes method. A particular feature of her approach 
is the formulation of the steady-state equations. It is effectively a node-based scheme, 
but in each cell the momentum and energy equations are assigned to the downstream 
nodes, and the continuity equation to the upstream nodes. This gives a stable scheme 
without the addition of any numerical smoothing, and it produces exceptionally good 
answers on very coarse grids. (It has recently been brought to the autlmrs’s attention 
that Wornam [3] has developed an implicit scheme for the quasi-one- dimensional Euler 
equations based on a similar idea. In Wornam’s scheme, for subsonic flow, the continu- 
ity and energy equations are included from the downwind direction and the momentum 
from upwind, and no numerical smoothing needs to be added.) 

Unlike most other upwdnd schemes that solve a Rieniann problem in order to obtain 
the numerical fluxes, the current approach computes the flux residual for a cell using a 
cell- vertex trapezoidal rule integration in exactly the same way as the cell- vertex schemes 
of Ni and Jameson (4]. The present scheme differs from those algorithms in that the 
distribution of the change in the state vector from the cells to the nodes is performed 
using a fully upwinded bias. The directions for the upwinding are based on the local flow 
direction, rather than being oriented with cell face normals, so that the upwinding is 
grid independent. Unlike other attempts at grid-independent upwind-differencing that 
require the construction of an upwind-biased difference stencil (e. g. Davis [5], Levy et 
al. [6]), the current scheme requires only information at the cell vertices to compute 
the distribution formulae for the cell. This compactness allows a simple data structure, 
with cell-to-node pointers being the only connectivity information required. Steady- 
state solutions are second-order accurate; however, for unsteady problems the scheme 
is only first-order accurate. 

Although the resulting scheme works quite well in predicting the pressure field, the 
lack of dissipation in the crossflow direction tends to decouple adjacent streamtubes, 
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resulting in odd-even oscillations in the density across streamtuhes. Also, large stag- 
nation density errors were observed near stagnation points. Two additional artificial 
viscosity terms were added to the scheme that smooth variations in entropy and vor- 
ticity. Although these viscosities are formally only first-order accurate, they actually 
improve the accuracy of the basic scheme in regions of the flow upstream of shocks, 
as they smooth quantities that should vanish analytically. This is at the cost of slight 
smearing of shock wakes. 

The scheme is described in the next section for a system of time dependent Unear 
hyperbolic partial differential equations in one space dimension. In Section 3, the exten- 
sion of the algorithm to the two-dimensional Euler equations is presented. An analysis 
of the steady-state difference operator for scalar advection is found in Section 4. The 
entropy and vorticity smoothing are explained in Section 5. Section fi presents airfoil 
results in two dimensions on a triangular mesh. Concluding remarks are made in Section 
7. 


2 One-Dimensional Scheme 


The standard first-order flux-upwinding algorithm for the linear system, 

ov du n 

"at 1 A dx r °’ 

is a cell-centered scheme. For a uniform mesh the semi-discrete form is 

r- (A ’ Uj + A" U 3+ \ ) - (A 4 Uj_i -i- A" Uj) = 0 

at 

where as usual 


with 


A TAT \ 

A :: 

diag(Xi), 
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. A + 

= diag(Xf ), 

A -TAT 1 

, A* 

diag( A t ~ ), 

A' max(0,A), 

A 

niin(0,A) 


This scheme can be rewritten using the following definition: 

sgn(A) = Tsgn(A)T~\ sgn(A) = dmj(sgn(Aj)) 


with 


sgn(A) -- < 


1 , A > 0 

1 , A < 0 

0 , A = 0. 


( 1 ) 

(2) 


Now then, 


(I f sgn(A)) A 
(I - sgn(A)) A 


2A ' , 
2A-, 


and hence 


A + - A(I f sgn(A))A, 

A = |(I - sgn(A))A. 

The difference equation may then be written 

+ 2 ft 1 + s S n ( A )) A ( u j - Uj ,) f (I sgn(A))A(U j+1 - Uj)] = 0. (3) 


Notice here the strong similarity to NTs distribution formulae. For a problem in 
which all characteristic speeds are identical, and using forward Euler time integration 
with a CFL number of 1.0, this scheme is identical to Nrs, which becomes purely 
upwinded under these conditions. 


Next, for the quasi-one dimensional Euler equations, 


dV 

dt 


, dF 


(4) 


the UCV (upwind control volume) scheme for a nonuniform mesh is 
( x j+i “ *j-i ((! » ■« n ( A i ..l))( F j F i-i - 

4 C 1 “ s f?n( A l))( F i+i - F i - (*j+ 1 - 
J+ 2 


)H, i) 

3 " 2 

,h m>) 


• (5) 


Notice that for irregular grids, this is an improvement over the usual flux-upwinding 
scheme that is only first order because now the source term H is being included with 
1 he correct value and the correct Ax so that in the steady state the solution is second 
order. Also, note that the scheme is conservative. This scheme for the one-dimensional 
Euler equations turns out to be identical to that proposed independently by van Leer, 
et al. [7]. 


3 Two-Dimensional Euler Scheme 


Now consider the full Euler equations of gas dynamics. In conservation form, the equa- 
tions are 

dU flF(U) <9G(U) 
dt ^ dx ^ dy 
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where 
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Here, p is the density; u and v y the velocity components in the x and y directions, 
respectively; p, the pressure; and E , the specific total energy. I he working fluid is 
taken to be an ideal gas with a constant specific heats ratio 7 , which gives 


In the semi -discrete 


p , _J_ P 
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Euler equations on an unstructured 


( 7 ) 


triangular mesh, the change in 


the state vector at a node has contributions from all of the cells which surround that 


Vi-%- - E(R,)k> <«> 

dt cell, 

\ j is the volume of the cell associated with node j, which can be defined as one third 
of the sum of the areas of the triangles around j, and (Rj)ic is the contribution to the 
flux residual at node j due to cell h . 

For a cell K with nodes 1,2,3 (numbered counter-clockwise) the flux residual Rr is 
defined in the usual node-based way: 


R* - ]rF Ay G Ar 

- \ ( F| (j/2 " 2/3) -Gi(*2-*3) 

f F 2 (»3-»l) - G 2(^3 * 1 ) 

4 F 3 (i/i y 2 ) G 3 (a:i - 12 ) )• ( 9 ) 


This residual is then distributed to the three nodes as follows: 

Ri«- - 5(1 - <r,Ani sgn(A,) + <r„AiiSgn(A„))RK, 

R 2 k j (I (r*An 2 sgn(A,) I <x„ A.s 2 sgn( A„)) Rr , (10) 

R 3K : | (I - (T,An 3 sgn(A„) I <r„Ai 3 sgn(A„)) Rr . 

a, and a n are parameters controlling the amount of upwinding in the streamwise and 
normal direction. Their usual values will be 1.0, while values of 0.0 will give Jameson s 
standard node-based algorithm. The streamwise and normal flux Jacobians A, and A n 
are defined by 

A. a - : A..?^ T Sy , 

A n A 11 X 1 Rrly, 


( 11 ) 


where A — 5F/9U and B — dG/dXJ. The unit, normal vector n — (n x ,n y ) r is defined 
to be normal to the unit streamwise vector s • (s T ,s y ) r which in turn is defined as 


_ tti + u 2 + u 3 

5 ~ .. ■ - ■ ■ . - 

|uj + U 2 -f U 3 | 

Finally, the terms A hj and A 5j are defined by these equations: 

An! = (£ 3 £ 2 ) ■ n, 

A n 2 ~ (£ 1 -x 3 )n, 

An 3 - (£ 2 *i)n, 


Asi = (£3 - £ 2 ) ■ s, 

As 2 = (£1 -£ 3 ) • 3, 

As 3 - (£2-£l)i*, 

An Anj 

Uj ' i(|An,| + !An 2 | + |An 3 |)’ 

A5 - 

" 2 (|Aji| -f |As 2 | + |As 3 |) (13) 

Thus the An/s and Aij’s vary from -1 to +1. 

4 Analysis of Two-Dimensional Scalar Algorithm 

To gain some understanding of the workings of the two-dimensional algorithm, consider 
its application to the scalar convection equation, 

du du 

HTT, "■ < 14 > 

which describes a wave moving in the positive x direction. To understand the behavior 
of steady-state solutions, consider the case when the time-derivative vanishes, i.e. 


First let us examine the fluxes distributed from the three equilateral triangles shown in 
Fig. 1. 
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The downwind nodes receive most of the flux residual, and the upwind nodes relatively 
little. Combining these results yields the steady- state flux operators for the meshes 
shown in Fig. 2. 


For grid 1, the steady-state equation is 
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For grid 2, the steady-state equation is 
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( 20 ) 


The best way in which to determine whether these discretizations are diffusive is to 
perform a Fourier analysis of the steady solutions. Define 

Ujk — z * exp (ik9) 


where 0 is the Fourier mode in the y-direction and z is the spatial growth factor in 
the t- direction. Because the grid is made of equilateral triangles, u jk must be carefully 
interpreted. The node numbering system is defined such that for even values of j the 
nodes are at integer values of fc, while for odd values of j the nodes are at the half-integer 
positions. This indexing system is indicated in Fig. 3. 

Substituting the Fourier mode into the steady-state equation for grid 1 yields a cubic 
equation for z: 

(z 2 l)(cos(0/2)z t 2) 0. (21) 
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The roots of this are 


2 


z = 1,-1, 


( 22 ) 


cos(0/2) 

The first two roots correspond to perfect convection of the quantity u along grid lines. 
This is possible because the convection velocity is perfectly aligned with the one set of 
edges of the computational domain. The third root is negative and has a magnitude 
which is greater than unity; this corresponds to an oscillatory, exponential type of 
solution. It is related to the fact that the discrete stencil includes nodes downstream of 
the central node, and it shows how the incorrect specification of the outflow boundary 
condition can lead to a local error which decays exponentially as one moves into the 
domain away from the boundary. 


Substituting the Fourier mode into the steady-state equation for grid 2 yields a 
quadratic equation for z: 

cos(0/2)z 2 + 2 (l + cos z (0/2)) z - 5cos(0/2) = 0. (23) 

One root of this quadratic is 

x _ -1 - cos 2 (fl/2) + y/(l + cpj 2 (0/2)) 2 +5cos 2 (fl/2) 

1 cos(0/2) 


5 cos(fl/2) 

1 + cos z (0/2) + ^(1 + cos z (0/2)) 2 4- 5cos 2 (0/2) 

5 


iW2) f cos ( 0 / 2 ) + 3 \/ 1 + 9*coj?fi/2) 


7 \ 2 r ( 2 ^) 

2 + - v^W2)) +3^1 + ^® 

From this last result it is clear that z x is real and lies in the range 0< z x < 1. Also, when 

6 < 1 , 

'> K TT¥m a 1 ■ 9</192 ' < 25 > 

This shows that the effect of the upwinding is comparable in nature to a fourth- difference 
smoothing when the computational grid is not perfectly aligned with the convection 
direction. 


The other root of the quadratic is 


z 2 = 


_5 

z i * 


(26) 


This root lies in the range z < -5 and so it again gives the oscillatory boundary-layer 
behavior at the outflow boundary. In an Euler equation application this would also give 
the exponential decay on either side of a shock. 
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5 Entropy and Vorticity Smoothing 


When applied to airfoil flows, the scheme as described above accurately predicts pres- 
sure distributions. However, two types of errors have been observed. The first is a 
significant error (~ 30 to 50%) in the stagnation density. It is highly localized, occur- 
ring at the grid point on the surface nearest the stagnation point. This error is not 
surprising, as the distribution formulae depend on the flow r direction which is singular 
at the stagnation point. The second error is due to the low dissipation in the crossflow 
direction, which tends to decouple adjacent streamtubes. This decoupling results in an 
odd-even oscillation of the density and velocity in the crossflow direction. The pressure 
field in any case is quite accurate and non- oscillatory, and no instabilities are observed. 

To fix these errors, two additional forms of artificial viscosity are added to the basic 
scheme: the first is a smoothing based on a second difference of entropy; the second 
is a smoothing based on vorticity. Both forms of the smoothing are conservative, and 
they have the property that they will return zero in an isentropic, irrotational flow. For 
airfoils in a uniform freestream, the entropy and vorticity smoothing do not corrupt 
the flow upstream of any shocks, and smear post- shock wakes only slightly. Also it has 
been found that the smoothing does not substantially smear the shocks. These artificial 
dissipation models were developed by the first author for his turbomachinery Euler code 
UNSFLO and are described in reference [ 8 ]. 

The entropy smoothing is based on a second difference in entropy. The (nondimen- 
sional) entropy for an ideal gas is defined as 


S = In — 7 In — , 

Poo Poo 

where the subscript oo refers to freestream conditions. Since it is desired that only 
those entropy variations that are decoupled from the pressure and velocity field are to 
be smoothed, the smoothing flux is based on 


5U e 


<9U 


dS 


p.« 


ss 


where d\J/dS is evaluated at constant pressure and velocity. This gives 


3U 

dS 



1 \ 

_ p 

u 

p,a 1 

V 


{uLtii) 


(27) 


Note that this is merely the projection of the right eigenvector of the flux Jacobians 
dF/d\J and dG/dJJ corresponding to entropy waves onto the state vector of con- 
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served quantities. This is simply a restatement that only entropy variations are being 
smoothed. 


The smoothing flux is then computed by taking the average of the entropy of the 
three nodes making up a cell, and distributing changes conservatively from each cell 
back to the nodes. The contribution of the smoothing flux from cell K to its t th node 
is then 


K 



{S K - so . 


(28) 


Here, (At/A)i is the timestep over the area for node i, (A/At)* is the area over the 
local timestep of cell if, 5* = (S\ + S 2 + S$)/3 is the average entropy over the cell, and 
t/ € is an arbitrary smoothing coefficient. 


Although this formula is formally only first order, in fact it improves the accuracy 
of the basic scheme upstream of shocks because the flow should be isentropic in that 
region. Thus any entropy variations in that region are due to the truncation error of the 
scheme. In shock wakes, of course, there exist entropy gradients, and the addition of 
the smoothing does smear the wake. For moderate values of v e this smearing is small. 

Inspection of Equation 27 shows that there is no smoothing of the crossflow momen- 
tum component. As a result, there will still be a decoupling of adjacent streamtubes, 
and a crossflow odd-even mode is allowed. To get rid of this error, a smoothing based 
on the vorticity is used. This smoothing is based on the observation that 


V x cJ = - V 2 u-bV(V*u) 


where w ~ V x u is the vorticity. Thus, the addition to the momentum equations of a 
term proportional to V x u? should smooth the flow when there are vorticity variations, 
but not corrupt irrotational flow. 


In two-dimensional flow, uj 


ku>, and the average vorticity in cell K is 



dK 


and the average value of V X 


is 


V x 


M - \ / 

dK 


a >ds. 


Based upon these observations, the vorticity smoothing is achieved by adding to the 
flux residuals of each node a term 


6 U? 


r 


j>&Fdy ^dx 

t dA l 


(29) 
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where AF and AG for the cell are defined as 

/ o \ 

v v pa curl (u) 


AF = - 


/ 


V 


o 

0 


AG = 


/ 


0 
0 

curl (u) 

0 


\ 


(30) 


where v v is an arbitrary smoothing coefficient, a is the speed of sound for the cell, and 
curl (u) is a scaled flow circulation, 

Au 2 iAx 3 i - Au 3 i A a;2i + At) 2 iAy 3 i - Au 3 iAy 2 i 
Y/AX21AJ/S1 - Ax 3 iAy 2 i 


curl (u) 


(31) 


where A () tj = {)i - {)j- The contributions of the vorticity smoothing fluxes of cell K 
to its surrounding nodes are 


S\] V lK = - y—J (AF K (j/3 - 2/2 ) AG k (*3 - * 2 )) , 

5U 2 K = ~ (~j~) (aFk (yi - y 3 ) AGx (x x - x 3 )) , ( 32 ) 

= - (^^ (AF* {y 2 -- J/i ) - AGk (x 2 - Xj)). 


The most attractive features of the entropy and vorticity smoothing are that they 
have a compact difference stencil, and they are well suited to transonic flows with 
a uniform freestream. However, these smoothing terms are not suited to high Mach 
number flows, in which there are strong shocks with large gradients of entropy and 
vorticity downstream of the shocks. For such flows, the smoothing errors are truly first 
order in most of the flow field. Also, the present scheme should be extended to the 
Navier- Stokes equations. The use of such smoothing for viscous flows would result in 
first-order errors thin-shear-layer regions. This latter problem may be moot, in that 
physical viscosity should suppress the odd-even mode and the stagnation density error 
as long as the viscous layers are properly resolved. In this case, there may be no need 
for the entropy and vorticity smoothing in the viscous regions of the flow. This question 
wilt be addressed in future work. However, the unsuitability of the present smoothing 
for inviscid flows with strong shocks remains. 

For flows with freestream Mach numbers significantly greater than 1, a background 
fourth- difference artificial viscosity similar to that of Jameson et al. [4] is preferable. 
Preliminary computations using fourth- difference smoothing have been run for airfoil 
flows from transonic to Mach 3.5 freestreams, and it works quite well; for transonic flows, 
the results are very similar to those obtained with the entropy and vorticity smoothing. 
Because the latter smoothing requires less memory, it is preferred for transonic flows. 
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6 Results 


To illustrate the capability of the scheme, solutions have been obtained for a standard 
inviscid test case given in [9], and referred to as AGARD 01. This case consists of 
flow over the NACA 0012 airfoil with a freestream Mach number M <*, = 0.8 and angle 
of attack a = 1.25°. The first set of results were obtained on the grid illustrated in 
Fig. 4. This grid was generated from a structured O-grid of 5248 nodes (128 X 41) 
by dividing each quadrilateral cell across a diagonal. The structured grid itself was 
used to obtain a benchmark solution using the proven code CFL2D (10), which uses an 
upwind-differencing algorithm. 

Figure 5 presents the surface c p distribution on the airfoil surface for CFL2D, the 
basic UCV scheme, and the UCV scheme with entropy and vorticity smoothing. Both <r, 
and <7 n were equal to 1, and the smoothing coefficients v e and v v for the third case were 
both taken to be 0.01. All three solutions agree very well with each other. The UCV 
scheme gives very sharp shocks, with only moderate pre- and post-shock overshoots. 
There are no undamped oscillations around the shocks, as will typically occur with 
central-differencing algorithms. Note that this is true of the solution without entropy 
and vorticity smoothing; the basic scheme UCV scheme gives nonoscillatory pressures. 
Also, note that the entropy and vorticity smoothing do not smear the shock significantly. 

From these results, it is not apparent why the entropy and vorticity smoothing are 
needed. To show why they are necessary, contours of constant density for the basic 
UCV scheme are presented in Fig. 6. Note the boundary-layer behavior at the wall, as 
well as the significant odd-even oscillations in the crossflow direction. These oscillations 
are particularly severe behind the shock, although they do not result in any instability. 
More serious, but not apparent in Fig. 6, is a large stagnation density error at the 
leading edge. The density at that point is 1.786, compared to the exact value of 1.351. 
Although not shown, the pressure contours are smooth, and the error in the stagnation 
pressure is much less severe, which indicates that the error lies in the failure of the 
scheme to compute the entropy correctly. 

With the addition of the entropy and vorticity smoothing (Fig. 7) both the boundary- 
layer and the odd-even mode are eliminated. The stagnation density now comes down 
to a more reasonable value, 1.388. It was found that the entropy smoothing alone failed 
to eliminate the odd-even mode, as there is no smoothing contribution to the crossflow 
momentum equation where the normal velocity vanishes. Both entropy and vorticity 
smoothing were found to be necessary. 

This test case has also been run using a background fourth- difference smoothing like 
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that of Jameson et al. [4] in place of the entropy and vorticity smoothing. A smoothing 
coefficient of 0.005 was found to be adequate to eliminate most of the stagnation density 
error (the computed stagnation density is 1.349 for this case). The surface c p distribution 
is shown in Fig. 8. Note that the primary difference between this result and the result 
obtained with entropy and vorticity smoothing is more smearing of the shocks, especially 
the very weak shock on the lower surface. Density contours are plotted in Fig. 9; they 
show that the fourth- difference smoothing also works well to eliminate the crossflow 
odd-even mode. 

The major reason for using an unstructured grid is the ability to handle arbitrary 
geometries. For complicated objects, it may be difficult to get a smooth grid. To 
illustrate how the scheme is suitable for distorted grids, a solution for AGARD 01 on a 
grid generated by Tim Barth and Dennis Jespersen (11] of NASA Ames Research Center, 
shown in Fig. 10, has been obtained. This grid has 6691 nodes, and was constructed 
by generating weighted pairs of pseudo-random numbers and Delaunay triangulation. 
Barth and Jespersen obtained a solution for the same test case on this grid, using their 
upwind algorithm presented in [11]. The UCV solution on this grid is shown in Figs. 11 
and 12. This solution was obtained with the entropy and vorticity coefficients set at 
0.01, as before. Again, the surface pressures are in good agreement with CFL2D. (Note: 
the CFL2D results are those obtained on the grid shown in Fig w 4, which has roughly 2.5 
times fewer nodes on the surface as the Barth grid.) Sharp shocks with small overshoot 
are obtained. The density contours in Fig. 12 are also seen to be very smooth, with no 
boundary-layer behavior. 


7 Conclusions 

A new Upwind Control Volume (UCV) scheme has been developed for obtaining numer- 
ical solutions to the Euler equations. This scheme has several very attractive features 
compared to existing methods. It is well suited to unstructured grids, providing ge- 
ometric flexibility. It is much simpler than most upwind schemes, in that it does not 
require complicated flux-limiting. The data structure required to implement the scheme 
on a triangular mesh is minimal: only cell-to-node pointers are needed. Although it is 
necessary to add smoothing to the basic UCV scheme in order to reduce stagnation 
density errors, the choice of entropy and vorticity smoothing does not corrupt the flow 
in irrotational, isentropic regions. The ability of the scheme to handle distorted grids 
without loss of accuracy has also been demonstrated. 
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Figure 6: Constant density contours, basic UCV scheme 



Figure 7: Constant density contours, entropy and vorticity smoothing 
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Figure 8: Surface Cp comparison, different smoothing 



Figure 9: Constant density contours, fourth- difference smoothing 
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Figure 12: Constant density contours, irregular grid 
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